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Tailored mixing inside individual droplets could be useful to ensure that reactions within micro- 
scopic discrete fluid volumes, which are used as microreactors in "digital microfluidic" applications, 
take place in a controlled fashion. In this article we consider a translating spherical liquid drop to 
which we impose a time periodic rigid-body rotation. Such a rotation not only induces mixing via 
chaotic advection, which operates through the stretching and folding of material lines, but also offers 
the possibility of tuning the mixing by controlling the location and size of the mixing region. Tuned 
mixing is achieved by judiciously adjusting the amplitude and frequency of the rotation, which are 
determined by using a resonance condition and following the evolution of adiabatic invariants. As 
the size of the mixing region is increased, complete mixing within the drop is obtained. 



PACS numbers: 47.51. +a, 47.61. Ne, 47.52.+j 



I. INTRODUCTION 



Droplets have been proposed as an alternative to stan- 
dard fluid-stream microfluidics for lab-on-a-chip appli- 
cations. This microfluidics approach, also referred to 
as "digital" because it uses "discrete" fluid volumes 
(droplets) rather than continuous streams, holds great 
promise due to the possibility of using single droplets as 
microreactors Efficient mixing, however, is needed 
for reactions to occur, but remains difficult to achieve 
because the Reynolds number (Re) is usually very small 
and so the flow is laminar. This issue has recently at- 
tracted much attention in the literature. For flows in 
microchannels, while there are many strategies based on 
altering the channel geometry, the use of forcing alone 
(see, e.g., 0, H, B H, M) has alsoproved to be efficient, 
especially in the case of low Re [7[ . The combination of 
both geometry alteration and forcing has been explored 
as well @, H, H, H3|- F° r droplet-based microfluidics, 
the forcing alone is the preferred strategy as the defor- 
mation of the droplet is difficult to control. In almost 
all cases, the enhancement of mixing in miniature ge- 
ometries is based on chaotic advection, the stirring phe- 
nomenon that stretches and folds fluid elements thus in- 
creasing the interfacial area between the two fluids to be 
mixed. Chaotic advection inside a liquid drop subjected 
to a forcing (at low Re) has been studied extensively 
[III [I2I . [H, Hj, [H, [H, H?} and was obtained experimen- 
tally by means of oscillatory flows [H, [I||. In this letter, 
we focus on unsteady - yet periodic - forcing. 
From a dynamical systems viewpoint, the introduction 
of a time-dependent perturbation or forcing breaks the 
invariants (related to the symmetries of the unperturbed 
system), thus introducing resonances between the nat- 
ural frequencies of the unperturbed problem and the 
frequency(-ies) of the forcing. Although such resonances 
create chaotic regions where mixing occurs, in general, 



chaotic and regular regions co-exist and unexpected reg- 
ular sizable pockets persist. 

In many situations where it is indeed possible to create 
chaos, controlling the mixing region(s) remains a chal- 
lenge. Such a control, however, should be possible since a 
chaotic system is sensitive to changes in parameter values 
(as it is to changes in initial conditions). These changes 
should generically modify the resonances, and thus the 
location and size of the chaotic regions. 
Our general approach along these lines is to consider a 
bounded three-dimensional (3D) flow, which is the super- 
position of an integrable flow vo with at least one invari- 
ant and a small time-dependent perturbation evi (x, t) , 
0<£<Cl. If Vo has only one invariant, the phase space 
contains two-dimensional tori. In this case, the perturbed 
flow, x. = vo(x) +evi(x, t), has poor mixing properties if 
the amplitude of the perturbation e is small, since two- 
dimensional (2D) tori act as barriers to chaotic diffusion 
(e.g., f20|). If, on the other hand, v has two invari- 
ants, trajectories of this integrable flow are all periodic. 
Most of these periodic orbits are expected to be broken 
by a generic perturbation vi with an arbitrarily small 
amplitude e. Efficient mixing properties might then be 
obtained with such perturbed flows. In this work, we 
consider an axisymmetric integrable flow possessing two 
invariants, thus possibly offering efficient mixing proper- 
ties after being perturbed. 

While many previous works [Il|, [IH, [2l], [22j have shown 
the existence of chaotic behavior in 3D bounded steady 
flows, we turn our attention to unsteady flows; the added 
unsteadiness targets the control of the chaotic behavior 
through resonance phenomena [I?], [H, [24| . Specifically, 
we seek to create a mixing zone of tunable size which re- 
mains localized within a well-defined region of the drop. 
This should also provide a rationale for the route to com- 
plete mixing as the perturbation increases. 
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II. 



MODEL 



A. Flow equation and assumptions 

We consider a spherical Newtonian drop immersed in 
an incompressible Newtonian flow in the case where the 
linear external field is characterized by translational ve- 
locity and vorticity vectors, similarly to [T3|. As in the 
latter reference, we assume that the local Re is much 
smaller than one and that the interfacial tension is suffi- 
ciently large for the drop to remain spherical. 
The internal velocity field is obtained by solving the 
Stokes flow problem for both the internal and external 
flows satisfying the continuity of velocity and tangential 
stress conditions across the drop surface. In addition, 
we introduce unsteadiness in the problem by making the 
vorticity time dependent. In a Cartesian coordinate sys- 
tem translating with the center-of-mass velocity of the 
drop, and with the z axis in the direction of the transla- 
tion, the paths of passive marker particles are given by 
the solution of the non-autonomous dynamical system: 



u = x = zx — a(t)uj z y, 

v = y = zy + a(t) (uj z x - lu x z) , 

w = z = (l-2x 2 - 2y 2 - z 2 ) + a(t)ui x y, 



(1) 



where all lengths and velocities have been non- 
dimensionalized by the drop radius and the magnitude 
of the translational velocity. Here, the vorticity is de- 
fined by u) = (u) x , u>y, uj z ) = (l/v2, 0, l/y/2), the uni- 
tary vector corresponding to the axis of rotation, and 
a(t) = e/2 (1 + cos ut), characterized by the frequency 
lj and the amplitude e. In this letter, we consider only 
small amplitudes, i.e. for < e <C 1. Note that the for- 
mer equations are identical to those in [l2T | except that 
the constant vorticity vector has been replaced by a(i)u). 
This can be done by either assuming unsteady vorticity 
in the external flow field, or by applying a time depen- 
dent body force. In practice, this could be realized, e.g., 
by creating a time dependent swirl motion in the external 
flow or by applying an electric field that exerts a torque 
on the drop (e.g.,[25[ it or work on electrorotation) . This 
flow is the superposition of a Hill's vortex and an un- 
steady rigid body rotation, and the surface of the drop, 



x 2 + y 2 



= 1, is invariant under flow {!]). 



B. Integrable case 



We now discuss some features of the unperturbed ax- 
isymmetric (2D) flow (e = 0). The flow possesses two 
independent integrals of motion, e.g., the streamfunction 
ip and the azimuthal angle <p: 



4' = l/2p 2 (l - r 2 ) , (/>= arctany/x, 



(2) 



where p 2 



y 2 and ip € [0,1/8]. The streamlines 
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FIG. 1: Streamlines inside the drop (without rotation) and 
their frequencies Q (ip) as given by Eq. ©. 



denoted by 1^ ^, and defined as (1 — 2p 2 ) 2 + (2pz) 2 = 1 — 
8ip (see Fig.[I|. Almost all streamlines are closed curves 
surrounding a circle of degenerate elliptic fixed points 
(p = 1/V2, z = 0). In addition, there are two hyperbolic 
fixed points located at the poles of the sphere which are 
connected by heteroclinic orbits. The frequency of the 
motion on is given by 
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(3) 

where j(ip) = \/T~~8tp and K is the complete elliptic 
function of the first kind. The frequency fl is bounded 
by two limits, fi(0) = and 0(1/8) = y/2 (see Fig. [[]). 

On every streamline T^,^, we introduce a uniform 
phase x rnod(27r) such that \ = on the x — y plane 
(with p < 1/V2) and \ = Q (ip). The unperturbed sys- 
tem, which can be rewritten in terms of (ip, 4>, \) as 

V^ = 0, = 0, x = Q(V), 
belongs to the class of action-action-anglc flows. 



C. Perturbed case 

In the perturbed case < e < 1, the time evolution of 
the two invariants of the unperturbed system is given by 



ip = -2a(t)oj x ip sin (f>G (ip, \) , 

4> = a(t)uj z — a(t)u> x cos 4>G (-0, x) 1 



(4) 



of the unperturbed system are lines of constant i/j and 



where G(V ) , x) = z/pis 2ir periodic in x an d has zero 
average in x- The time evolution equation for x is 

X = n(4,)+a(t)H(ij,0,x), 

where H is 2ir periodic in x- The dynamics possesses two 
time scales, a fast one (of order one) associated with x, 
and a slow one (of order 1/e) associated with ip and <p. 
If f2 and to are incommensurate, then the averaging over 
51 and over to can be performed independently. In this 



case, the time-periodic terms in Eq. (J4j) average out, and 
the averaged system reduces to ip = 0, <j> = —e/2. Thus 
in the averaged system the value of ip is conserved as it 
was in the unperturbed system; in other words, ip is an 
invariant of the averaged system. Each trajectory of the 
averaged system evolves on two-dimensional nested tori 
In the perturbed system, ip is an adiabatic invariant 
and the motion follows adiabatically the tori T^. 

III. METHODS AND RESULTS 

A. Mixing generation via resonance phenomena 

We now turn to the generation of a 3D chaotic mixing 
region inside the drop, for which we seek to control both 
the location and the size. The strategy used for this 
purpose is to bring a chosen family of unperturbed tori 
7^, into resonance with the perturbation a(t) by adjusting 
the frequency lo to satisfy the resonance condition 

nfl(ip) - lo = 0, (5) 

for some n £ N (see Fig. [TJ. For any fixed lo we denote 
by {T^ n \uj) \n <G N} the set of resonant tori %p satisfy- 
ing ([5])- Hereafter, we denote the chaotic mixing region 
generated around T^ x \lo) by CMR. 

B. Control of the mixing 

FiguresOand Opresent Liouvillian sections of the per- 
turbed system, which consist of 2D projections of time- 
periodic 3D flows by a combination of a stroboscopic map 
and a Poincare section (here, the y = plane). Figure [2] 
shows that a perturbation a(t) creates a 3D CMR around 
T' 1 ) (lo) and its location is controlled by varying lo accord- 
ing to Eq. (O. In what follows, we analyze the location 
and the size of the CMR as to and e vary. 

For small values of w, all resonances are located near 
the pole-to-pole heteroclinic connections (at ip = 0, near 
the z axis and near the boundaries of the drop, see 
Fig. [2^i). As lo is increased, the CMR penetrates deeper 
into the drop (Fig. [2ji>) . In the interval < lo < \/2, the 
CMR is the largest chaotic region (compared to chaotic 
regions corresponding to higher order resonances), with 
all the other chaotic regions localized close to the z axis 
and near the drop boundaries (around the heteroclinic or- 
bits); this is due to the shape of f2 (ip). As lo is increased 
further, the CMR moves toward the location of the ellip- 
tic fixed points of the unperturbed system, closely follow- 
ing the location of the resonant torus T^(lo) (Fig. [2};). 
As the value of lo approaches the CMR shrinks to 
the circle of elliptic fixed points (Fig. [5]). 

Whereas the frequency lo of the rigid body rotation is 
mostly responsible for the location of the CMR, it is its 
amplitude e which mostly determines its size. Figure [3] 
shows that the size of the chaotic mixing regions created 




FIG. 2: Liouvillian sections for the amplitude e = 0.03 and 
the frequencies w = 0.55,0.93,1.28,1.41 (a-d). The (red) 
dashed line inside the CMR is the torus T^K 




FIG. 3: Liouvillian sections for the frequency lo — 1.376 and 
the amplitudes e = 0.01, 0.05, 0.10, 0.20 (a-d). 



by the n = 1 resonance and by higher order resonances 
(mostly the n = 2 resonance) increases as the ampli- 
tude of the perturbation increases. Around e 0.20, the 
chaotic regions around the heteroclinic orbits and the 
CMR join together to cover the entire drop volume. 

Recall that in the averaged system the adiabatic in- 
variant ip is constant. In the exact system, however, 
along a given trajectory starting at ip = ipo it varies be- 
tween ip~ (tp ;LO,e) and tp + (ipo]LO,e). The width Aip = 
ip + (ipo; lo, e) — ip~ (ipo;uj, e) is small away from the reso- 
nance, and increases significantly closer to the resonance. 
The projection of three characteristic trajectories onto 
the (ip, </>)-plane (called the slow plane in dynamical sys- 
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FIG. 4: Projection of three characteristic trajectories on the 
slow phase plane, with (pa — and ipo = 0.010, 0.073, 0.125. 
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mal CMR size one can reach for a given amplitude e of 
the rotation. On the other hand, Atp versus e increases 
quite monotonically for all values of lo. The derivation 
of the maxima locations and estimates of Atp as func- 
tions of the parameters and the order of resonance, will 
be addressed elsewhere. The structure of the CMR in 
our case is rather different from that obtained in other 
problems that possess resonance-induced chaotic advec- 
tion. Namely, here the size of the CMR vanishes as e 
goes to and the CMR is localized near the resonance. 
In contrast, in the flow considered in, e.g., fl7l |. the mix- 
ing is caused by resonances, but the CMR occupies a 
volume on the scale of the whole system. The difference 
comes from the fact that the averaged change of the fre- 
quency of the fast system vanishes in the current system, 
thus preventing the trajectories starting away from the 
resonance from approaching it. This property makes the 
kind of flows investigated here useful as it may be ad- 
vantageous to localize the mixing in certain parts of the 
system only. 



0.04 



0.08 



0.12 



0.16 E 0.20 



FIG. 5: Size of the chaotic mixing region; Upper panel: Nor- 
malized Aip vs. lo for the amplitudes e = 0.01, 0.05, 0.10, 0.20 
(a-d); Lower panel: Normalized Alp vs. e for u = 
0.55,0.93,1.28,1.41 (a-d). 

terns) is presented in Fig. 0] The narrow regions on the 
sides are off-resonance trajectories that stay quite close 
to the corresponding tori %p. In between, the middle 
trajectory deviates much further from its %/, = T^'iui) 
and fills the entire CMR. The quantity A?/> is probably 
the most convenient quantity to estimate the size of the 
CMR (around T^> (lo) for lo < The volume between 
the tori (lo) and T^+(lo) gives the CMR size in 3D. 
The dependence of the size of the CMR (in terms of Aip) 
on e and w is illustrated in Fig. [5] The curves (a)-(d) in 
the upper and lower panels correspond to the Figs. ^.-d 
and [3£i-d, respectively. For a given lo value (i.e. for a 
given T^(lo)), the size can be controlled by adjusting 
the value of e; for example in the range of frequencies 
1.181 < u> < 1.357, the entire droplet exhibits chaotic 
mixing for e > 0.175. For each smaller value of e the 
size reaches a maximum for a certain value lo" 1 (e) of 
the frequency. On the one hand, this property can be 
used as an optimization technique to obtain the maxi- 



IV. CONCLUSION 

In summary, we have shown that by applying a judi- 
cious oscillatory rotation to a translating drop (an inte- 
grable system), one can create a chaotic mixing zone with 
a prescribed location and size. The appropriate values of 
the parameters of the perturbation (here, a rotation of a 
given frequency and amplitude) are determined by quan- 
titative features of the integrable case. For any amplitude 
of the rotation, the frequency optimizing the CMR size 
has been obtained. Such an optimization could be useful 
in guiding the design of practical mixing devices aiming 
at the best possible mixing rate within individual drops. 
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